A minimizing principle for the Poisson-Boltzmann equation 
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The Poisson-Boltzmann equation is often presented via a variational formulation based on the 
electrostatic potential. However, the functional has the defect of being non-convex. It can not 
be used as a local minimization principle while coupled to other dynamic degrees of freedom. We 
formulate a convex dual functional which is numerically equivalent at its minimum and which is 
more suited to local optimization. 
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Introduction 

The Poisson-Boltzmann treatment of free ions is widely 
used in simulations with implicit solvents. It successfully 
describes the mean field screening properties of ionic so- 
lutions. At the most phenomenological level it is found 
by balancing two crucial features of the ionic system: The 
electrostatic energy coming from Coulomb's law, plus the 
ideal entropy of mixing of the ions. Written in terms 
of ion concentrations Cj , or the charge concentrations 
Pj = Cj q, we find the free energy of a solution from 
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where the total charge density p = Y.] Pj + Pf ■ Pf is an 
external fixed charge density - associated with surfaces 
or molecular sources, Cjo is the reference density of com- 
ponent j. If we minimize this functional of Cj then we 
find an effective free energy for the source field p f , [1]. 
The problem in this formulation is the appearance of the 
long-ranged Coulomb interaction [2] which renders the 
evaluation and minimization less efficient than might be 
wished. 

We start with a few points of notation: We will use 
the symbol f to describe a density of free energy of an 
electrolyte in mean field theory, h denotes a complex ac- 
tion which occurs in a functional integral and u is the 
electrostatic energy of the electric field. We will freely 
integrate by parts in our expressions dropping boundary 
terms. The transformations that we perform conserve 
the stationary point of the mean field solution which is 
denoted f s . The statistical weight of a configuration is 
then 
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no matter the arguments of the functional f . Pf is kept 
constant during all variations of the electric and ionic 
degrees of freedom. We denote the Legendre transform 
of a convex function g(x) as 



g(t)=£(g(x))[£J=sup(x£.-g(x)) 



(3) 



In eq. (1) one conventionally decouples the electrostatic 
interactions by introducing the potential as an additional 



variational field. If we do so then we find that 
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f = pcj> — e- 
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We see at once that we have gained in locality of the 
formulation; the other advantage is that it is valid for 
arbitrary spatial variation in the dielectric properties, 
e(r). Unfortunately, the counterpart is that the result- 
ing free energy is no-longer convex. We do not have a 
minimizing principle, rather a stationary principle. This 
excludes some of the simplest optimization strategies that 
one might want apply - such as simultaneous annealing 
of conformational and electrostatic degrees of freedom. 
From eq. (4) one studies the differential 

df = Y_ dcj (qj<p + k B Tln(cj/Cj )) + dcb(p + V • eVcb) 



We now impose that the coefficient of dCj is zero so that 



qjcf) + k B Tln(cj/Cj ) = 



(5) 



Substituting back into eq. (4) we find the standard form 
for the Poisson-Boltzmann functional [3-71 : 



f = p f cp - e 



k B T^Cj e q 



(6) 



We take the variations of this functional with respect to 
<t) to find: 



p f + V • eV4> + qjCj e 







For the symmetric ion system 
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f =p f 4> - e- 
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p f + V • eVdp - 2qc sinh(|3q4>) = 



(7) 



(8) 
(9) 



Eq. (8) continues to be awkward numerically since clearly 
both the derivative and the cosh functions are unbounded 
below; simple annealing procedures are thus unstable. 
One has to solve the partial differential equation eq. (9) 
in practical applications [8, 9]. 

The purpose of this paper is to derive functionals that 
are equivalent to those given above which combine the 



2 



advantages of the convexity of eq. (1) and the locality of 
eq. (4). We show how to find functions which are numeri- 
cally equivalent to those widely used in the literature. We 
consider this absolutely crucial since much time has al- 
ready been invested in calibration of potential functions. 
We thus look for ways of rendering convex known, ac- 
curate functionals. Our main tool in this effort will be 
the Legendre transform in a form presented in [10] as a 
reciprocal principle for variational calculations. In the 
following calculations the way that the energy e(Vcb) 2 /2 
is transformed into an equivalent form D 2 /2e will remind 
the reader of the transformation of the kinetic energy in 
a Lagrangian, rrvx 2 /2 into the Hamiltonian equivalent, 
p 2 /2m. This is a standard example of the use of Legen- 
dre transforms in classical mechanics. 

We start by summarizing the state of the art in the cal- 
culation of effective free energies for Poisson-Boltzmann 
free energies, and then showing in a simple example how 
the ideas of dynamically constrained fields can be used 
to write a convex free energy in terms of the electric dis- 
placement field D. We then generalize the procedure to 
arbitrary formulations of the free energy demonstrating 
the link to the theory of Legendre transforms. 



Field theory route to Poisson-Boltzmann 



Constrained fields plus mixing entropy 

An alternative, local approach to electrostatics was in- 
troduced in [16-18] to avoid the use of Ewald and fast 
Fourier methods in electrostatic codes. The ideas were 
applied to the Poisson-Boltzmann equation in [19, 20]. 
A convex energy is associated with the Gauss law con- 
straint which is implemented dynamically. In particular 
we reproduce the mean field Poisson-Boltzmann equa- 
tions from the functional 

f =^+k B T^(c j ln(c j /c j0 )-c j ) (11) 
j 

with divD-p = (12) 

In the case of a one-component plasma one trivially elim- 
inates the density degree of freedom to find 

D 2 

f = ^+k B Ts((divD-p f )/q) (13) 
ze 

s(z) = zln(z/co) — z (14) 

The functional eq. (13) is now an unconstrained func- 
tional of the field D. It can be discretized as in [16] 
where the vector field D is associated to the links of a 
cubic lattice, while scalar quantities such as (divD) and 
pf are associated with the vertexes of the lattice. Eq. (13) 
has the advantage of being both convex and local. 



In recent years theoretical methods have been intro- 
duced to generalize the application of Poisson-Boltzmann 
equations to a larger range of systems, [11, 12]. All 
these methods are based on the Hubbard-Stratonovich 
transformation to break pair interactions into one-body 
potentials. They have shown their power in producing 
functionals that include the finite volume of ions as well 
as mobile dipoles. The Hubbard-Stratonovich transfor- 
mation produces a complex action, which can however 
be studied at its saddle point. In this case one repro- 
duces the known Poisson-Boltzmann functionals for sim- 
ple situations, and in addition study fluctuations about 
the mean field solution in terms of a loop expansion. 

A typical example of such a function for both symmet- 
ric ions and dipoles is [13-15] 
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where At on and A dip are the chemical potentials of the 
ions and the dipoles, a a hard core size, po a dipole mo- 
ment. e(r) = e , since dielectric effects are generated 
dynamically. The logarithmic term in the energy is found 
from studying the partition function of a lattice gas. 



Reciprocity and Legendre transforms 

We now show how to transform a Poisson-Boltzmann 
functionals expressed in terms of electrostatic potentials 
cb into those based on the electric fields, D. Start with 
the electrostatic energy density 



u = pf qb — e 
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As is well known the Poisson equation is found as the sta- 
tionary condition of eq. (15). Introduce now E = —Vet) 
using a Lagrange multiplier D to find the constrained 
functional 



E 2 

u = p f 4> - e— + D • (E + Vcb) 



(16) 



Integrate by parts and regroup: 



E 2 



u = -e— +D-E-4>(divD-p f ) (17) 



We now eliminate E and see the equivalence to 



D 2 



u= 4>(divD - p f ) 

2e 



(18) 



We now re-interpret the field cb as a Lagrange multiplier 
[10] imposing Gauss' law and find the starting point of 
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our previous papers [16] on constrained statistical me- 
chanics in electrostatics. 

Let us now generalize this approach to Poisson- 
Boltzmann functionals expressed in terms of the poten- 
tial. 

f =p f cb-e^^-g(<j>) (19) 

We will explicitly treat the example of the symmetric 
electrolyte solution where g(4>) = 21cbTco cosh ((3q<J)). 
The transformation of eqs. (15-18) still goes through and 
we find 

f=5^+*(Pf-divD)-g((|>). (20) 

Further variations with respect to (f> give simply the Leg- 
endre transform, eq. (3), of the function g((b), where the 
transform variable is £, = (pf — div D) 

To continue the transformation we require the station- 
ary point of 

4>£,-2k B Tc cosh(|3qcb) (21) 
finding ct> = (k B T/q) sinh -1 (£,/2qco) and 

£(g)H = sinh- 1 (f,/2qc ) - k B Ty^ ~ ^ 2/q 2 

(22) 

For symmetric solutions we conclude that the reciprocal 
functional that we require is 

D 2 

f = — +£(g)[p f -divD] (23) 

This demonstrates the principle result of the present pa- 
per. Starting from a standard, concave functional ex- 
pressed in terms of the potential we have found a convex 
functional of the vector field D. The critical points of 
the two functionals are, however, numerically identical. 

When there are no ions within a region - for instance 
within a macromolecule - the Legendrc transform re- 
quires some care in its definition: We take g(<$>) — r\4> 2 /2 
with T| small. Then the Legendre transform is g = 
£, 2 /(2r|). Using finite r| gives a Yukawa or Debye in- 
teraction, with screening of the electrostatics. The limit 
r) — > imposes a delta-function constraint on Gauss' law. 
Experience with local Monte Carlo algorithms based on 
the electric field [16, 17] implies that relaxation of longi- 
tudinal and transverse degrees of freedom via link, and 
plaquette updates is the most efficient manner to sample 
the above functionals. 



Dipolar systems 

We now consider the transformation of the functional 
eq. (10) for small chemical potentials Aj , where we expand 



the logarithm to lowest order. In this case we find a 
contribution of the form 

E 2 sinh(ftpo|E|) 
f E — eo T - A d ppo|E| ( 24 ) 

Again if we consider the extremal equation for E we rec- 
ognize the Legendre transform, this time for the electric 
field, rather than for the potential. We have here defined 
A,- =k B TAj/a 3 . 

Thus the fully transformed functional in presence of 
both ions and dipoles is 

/ E 2 - anh(gpo|E|) \ 
f = ^e 0T+Adip ppo|E| j[D]+ 

2A ion £(cosh((3e4>))[divD-p f ] (25) 

For small fields, E, we expand the first line of eq. (25) 

and find eo^-(1 + A^ip (3 2 p 2 /3eo). The modification of 
the curvature of the function is a manifestation of the 
electric susceptibility the dipoles. 

Evaluating transforms 

In general it is impossible to analytically transform the 
functions needed in the dual or reciprocal formulation. In 
particular the full functional eq. (10) is not a simple sum 
of terms as assumed in the derivation of eq. (25). How- 
ever in a given ionic system the Legendrc transformed 
functionals arc uniform in space. They should be cal- 
culated just once for a given set of chemical potentials. 
For the general problem we require the two-dimensional 
transform with respect to the variables (cb, E). Standard, 
fast numerical methods [21] are available for performing 
these transforms, and the result can be tabulated for a 
given application. For the highest accuracy an iterative 
Newton-Raphson step can be used to converge the in- 
terpolated results to very high accuracy. Forces can be 
evaluated as well as energies, for instance in eq. (22), 
dg(£j/d£, = <fy - a general property of the transform [22], 
thus the force on a test particle with charge is just 

F ) =-e ) V<|) = ejE (26) 

Link with Infimal Convolution 

The Legendre transformation of the entropy function 
s(z) of eq. (14) is the Boltzmann factor 

£(s(z))H =c e £ (27) 

appearing in eq. (6). The fact that the effective action 
eq. (6) is a sum of exponential functions shows that the 
reciprocal action must be linked to infimal convolution of 
the original free energy defined as follows [23] : 

£(g + d)[x] = inf [g(y) + d(x— y)] (28) 
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We now show why this is the case for the symmetric 
electrolyte. 

Consider a symmetric electrolyte with unit charges, 
q = 1 and with reference concentrations Co = 1 • We im- 
pose Gauss' law with a Lagrange multiplier <J> and study 
the constrained minimum of the total ionic entropy 

t(d , c 2 ) = s(ci ) + s(c 2 ) + 4>(t - ci + c 2 ) (29) 

where f, = (divD — Pf). The variation with respect to Cj 
gives 



(30) 



We now study this Legendre transform of the sum of two 
functions. For the symmetric electrolyte we require (from 
the above definition of the infimal convolution) 

t'(« = inf [y ln(y) -y + (y - £,) ln(y - E.) - [y - £,)] 
v 

(31) 

giving y(y-£J =1,y = (£,+V4 + f, 2 )/2. Thus the action 
for the symmetric electrolyte eq. (22) is found as the in- 
fimal convolution of two ideal entropy functions with the 
help of the identity sinh -1 (£/2) = In ((£, + V4 + £, 2 )/2) 



The only thing that remains is the integral over cb. This 
we recognize as a Fourier transform with variable f, = 
(divD — pf ). Thus the action expressed in terms of D is 

D 2 

h= — -ln{J-(e-9 (tt>) )}[divD-p f ] (38) 

with T the Fourier transform. 

This action should be integrated over to find the par- 
tition sum 



Z= dDe-^ hd3r 



(39) 



When the fluctuations at the saddle point are neglected 
the Fourier transform reduces to the Legendre transform 
as above, and the infimal convolution is a simple echo of 
the standard convolution of statistical weights occurring 
at the saddle point. If the Fourier transform becomes 
negative then the action becomes complex, or one must 
at least sample functions which are non-positive. Again 
in general we do not expect to be able to evaluate the 
Fourier transform analytically but again tabulation for a 
specific problem is possible. 



Fluctuations and Fourier transforms 



Conclusions 



We now return to the full field theoretical formula- 
tion of the Poisson-Boltzmann equation after Hubbard- 
Stratonovich transformation of the action, now denoted 
h, but before the saddle point evaluation [11, 12]. As a 
concrete example we consider the symmetric electrolyte: 

h=e l -2k B Tc cos((3q4)) ip f cb (32) 



—e 



2 

(V4>) 2 



+ g(4>) -iPf<t> 



(33) 



We now consider transformation in the philosophy of 
our above reciprocal formulation, but replacing Lagrange 
multipliers by complex integral representations of the 
delta-function. We do not neglect fluctuations in the 
fields and do not make any approximation in the statisti- 
cal mechanics of the ionic system. Following very closely 
the logic described above we find the following succession 
of transformations of the action: 

H=e^^ + g(4))-ip f 4> (34) 
E 2 

-^e— + g(ct>)-iD-(V4> + E)-ip f 4> (35) 
E 2 

^e— + g(ct)) -iD • E + icJ)(divD - p f ) (36) 
D 2 

^27 + 9(4>)+i4>(divD-p f ) (37) 



To conclude, we have introduced a duality transfor- 
mation for Poisson-Boltzmann functionals which allows 
us to find a local minimizing principle for both electro- 
static and conformation degrees of freedom in a simu- 
lation. This opens the perspective of simpler annealing 
and dynamic relaxation in molecular simulation, includ- 
ing local Car-Parrinello evolution of ionic degrees of free- 
dom [17] for complex molecules. Clearly the interpola- 
tion of source charges to a grid requires control of their 
self-energy in a manner which in familiar in molecular 
dynamics simulations [24-26]. 

The functionals are designed to exactly conserve the 
saddle point. This is only true in practice if the numer- 
ical discretization is strictly equivalent. The final for- 
mulation requires a discretized divergence for Gauss' law 
in eq.(20), "Div". This divergence operator is adjoint 
to the "— Grad" from eq. (16). Finally the Laplacian in 
the potential formulation must obey V 2 = Div Grad. If 
this is not true discretization errors differ between the 
formulations. 

The original and reciprocal functionals have very dif- 
ferent forms- as noted above a weakly constraining en- 
ergy r|cp 2 /2, with r\ small, is transformed into a steep 
constraint £, 2 /(2r|), and vice- versa. In the future it will 
be interesting to understand how this influences mathe- 
matical descriptions of fluctuations about the mean field 
behaviour. As stated in [10] the potential and field for- 
mulations can be used together to give upper and lower 
bounds for the mean field free energy. 
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